# Timed version of the main run function with GPTL-style instrumentation

using .Timing: @timed_function, reset_timing!, print_timing_report

"""
    run_timed(settings_files::AbstractString, settings_config::AbstractString)

Run FLEXINVERT with detailed timing instrumentation similar to Fortran GPTL.
"""
function run_timed(settings_files::AbstractString, settings_config::AbstractString)
    reset_timing!()

    # Main timing block
    result = @timed_function "main" begin

        # Initialization
        files, config = @timed_function "initialization" begin
            f = read_file_settings(settings_files)
            c = read_config_settings(settings_config)
            (f, c)
        end

        # Initialize Domain
        domain = @timed_function "init_domain" begin
            Domain(config)
        end

        # Initialize transport system
        grids, date_info = @timed_function "init_transport" begin
            try
                initialize_transport_system(files, config, domain)
            catch err
                @warn "Transport initialization failed" err
                (nothing, nothing)
            end
        end

        # Read orography (placeholder - often skipped)
        @timed_function "read_orog" begin
            # Usually empty in our setup
            nothing
        end

        # Read land-sea mask
        @timed_function "read_lsm" begin
            if !isempty(files.file_lsm) && !isempty(files.varname_lsm)
                try
                    NetCDFIO.read_lsm!(domain, files.file_lsm, files.varname_lsm,
                                      files.lonname_lsm, files.latname_lsm)
                catch err
                    @warn "Failed to read LSM" err
                end
            end
        end

        # Get regions
        @timed_function "get_region" begin
            get_regions!(domain, files, config)
        end

        # Read receptors and observations
        receptors, obs = @timed_function "read_obs" begin
            recept = Vector{Receptor}()
            observations = empty_observations()

            if !isempty(files.file_recept) && isfile(files.file_recept)
                recept = read_reclist(files.file_recept)
            end

            if !isempty(files.path_obs) && isdir(files.path_obs) && !isempty(files.suffix)
                observations = read_observations(files.path_obs, files.suffix;
                                                datei=config.datei, datef=config.datef)
                if config.measerr > 0
                    apply_measurement_error_floor!(observations, config.measerr)
                end

                # Compute initial condition contributions (cini/cinipos)
                if length(observations.concentrations) > 0
                    @timed_function "init_cini" begin
                        try
                            compute_initial_contributions!(observations, files, config, domain)
                        catch err
                            @warn "Failed to compute initial condition contributions" err
                        end
                    end
                end
            end

            (recept, observations)
        end

        # CO2/GHG processing (includes covariance setup)
        prior_fluxes = nothing
        covariance = nothing
        footprints = nothing

        if length(obs.concentrations) > 0
            prior_fluxes, covariance = @timed_function "co2_ghg" begin
                # Read prior fluxes
                fluxes = read_prior_fluxes(files, config, domain)

                # Build covariance with sub-timings
                cov = @timed_function "error_cov" begin
                    build_prior_covariance(config, domain, fluxes)
                end

                # Correlation computation (if needed)
                @timed_function "correl" begin
                    # Correlation calculations are embedded in covariance build
                    nothing
                end

                # Perturbation (usually not used in analytic)
                @timed_function "perturb" begin
                    # Perturbation for ensemble methods
                    nothing
                end

                (fluxes, cov)
            end

            # Load footprints
            footprints = @timed_function "load_footprints" begin
                load_observation_footprints(files, config, domain, obs)
            end
        end

        # Average footprints (if enabled)
        @timed_function "average_fp" begin
            if config.average_fp && length(obs.concentrations) > 0
                try
                    average_footprints!(obs, config, files)
                catch err
                    @warn "Footprint averaging failed" err
                end
            end
        end

        # Retrieval/Inversion
        inversion_result = nothing
        state = nothing

        if length(obs.concentrations) > 0 && prior_fluxes !== nothing && covariance !== nothing
            state, inversion_result = @timed_function "retrieval" begin
                # Initialize state
                st = initialize_state(config, domain, covariance)

                # Run forward model (simulate)
                y_mod = @timed_function "simulate_analytic" begin
                    forward_model(config, domain, prior_fluxes, obs;
                                footprints=footprints, files=files)
                end

                # Inversion calculation
                result = @timed_function "analytic" begin
                    if config.method == "analytic"
                        invert!(st, prior_fluxes, obs, covariance, config, files, domain, footprints)
                    else
                        nothing
                    end
                end

                # Product of inverse covariance (prod_icov)
                @timed_function "prod_icov" begin
                    # This is typically called within invert!
                    nothing
                end

                @timed_function "prod_icov" begin
                    # Called twice in Fortran
                    nothing
                end

                (st, result)
            end
        end

        # Save files
        @timed_function "save_file" begin
            if !isempty(files.path_output) && inversion_result !== nothing
                try
                    save_optimization_logs(inversion_result, files)
                    OutputWriter.save_region_products(files.path_output, domain)
                catch err
                    @warn "Failed to save outputs" err
                end
            end
        end

        # Return results
        (; files, config, domain, receptors, obs, grids, date_info,
           prior_fluxes, covariance, state, inversion_result, footprints)
    end

    # Print timing report
    print_timing_report()

    return result
end

export run_timed
